The smuthi.linearsystem package

linearsystem

This package contains functionality that is related to the assembly or solution of the system of linear equations that yield the solution of the scattering problem.

linearsystem.linear_system

This package contains classes and functions to represent the system of linear equations that needs to be solved in order to solve the scattering problem, see section 3.7 of [Egel 2018 dissertation].

Symbolically, the linear system can be written like

\[(1 - TW)b = Ta,\]

where \(T\) is the transition matrices of the particles, \(W\) is the particle coupling matrix, \(b\) are the (unknown) coefficients of the scattered field in terms of an outgoing spherical wave expansion and \(a\) are the coefficients of the initial field in terms of a regular spherical wave expansion.

class smuthi.linearsystem.linear_system.CouplingMatrixExplicit(vacuum_wavelength, particle_list, layer_system, k_parallel='default', use_pvwf_coupling=False, pvwf_coupling_k_parallel=None)

Class for an explicit representation of the coupling matrix. Recommended for small particle numbers.

Parameters:
  • vacuum_wavelength (float) – Vacuum wavelength in length units
  • particle_list (list) – List of smuthi.particles.Particle objects
  • layer_system (smuthi.layers.LayerSystem) – Stratified medium
  • k_parallell (numpy.ndarray or str) – In-plane wavenumber. If ‘default’, use smuthi.fields.default_Sommerfeld_k_parallel_array
class smuthi.linearsystem.linear_system.CouplingMatrixPeriodicGridNumba(initial_field, particle_list, layer_system, periodicity, ewald_sum_separation_parameter, num_threads='default')
Class for an explicit representation of the coupling matrix of periodic particle arrangements.
Computation supports Numba.
Parameters:
  • initial_field (smuthi.initial_field.PlaneWave) – initial plane wave object
  • particle_list (list) – list of smuthi.particles.Particle objects
  • layer_system (smuthi.layers.LayerSystem) – stratified medium
  • periodicity (tuple) – (a1, a2) lattice vector 1 and 2 in carthesian coordinates
  • ewald_sum_separation_parameter (float) – Ewald sum separation parameter
  • num_threads (int or str) – if ‘default’ all available CPU cores are used if negative, all but num_threads are used
class smuthi.linearsystem.linear_system.CouplingMatrixRadialLookup(vacuum_wavelength, particle_list, layer_system, k_parallel='default', resolution=None)

Base class for radial lookup based coupling matrix either on CPU or on GPU (CUDA).

Parameters:
  • vacuum_wavelength (float) – vacuum wavelength in length units
  • particle_list (list) – list of sumthi.particles.Particle objects
  • layer_system (smuthi.layers.LayerSystem) – stratified medium
  • k_parallel (numpy.ndarray or str) – in-plane wavenumber. If ‘default’, use smuthi.fields.default_Sommerfeld_k_parallel_array
  • resolution (float or None) – spatial resolution of the lookup in the radial direction
class smuthi.linearsystem.linear_system.CouplingMatrixRadialLookupCPU(vacuum_wavelength, particle_list, layer_system, k_parallel='default', resolution=None, interpolator_kind='linear')

Class for radial lookup based coupling matrix running on CPU. This is used when no suitable GPU device is detected or when PyCuda is not installed.

Parameters:
  • vacuum_wavelength (float) – vacuum wavelength in length units
  • particle_list (list) – list of sumthi.particles.Particle objects
  • layer_system (smuthi.layers.LayerSystem) – stratified medium
  • k_parallel (numpy.ndarray or str) – in-plane wavenumber. If ‘default’, use smuthi.fields.default_Sommerfeld_k_parallel_array
  • resolution (float or None) – spatial resolution of the lookup in the radial direction
  • kind (str) – interpolation order, e.g. ‘linear’ or ‘cubic’
class smuthi.linearsystem.linear_system.CouplingMatrixRadialLookupCUDA(vacuum_wavelength, particle_list, layer_system, k_parallel='default', resolution=None, cuda_blocksize=None, interpolator_kind='linear')

Radial lookup based coupling matrix either on GPU (CUDA).

Parameters:
  • vacuum_wavelength (float) – vacuum wavelength in length units
  • particle_list (list) – list of sumthi.particles.Particle objects
  • layer_system (smuthi.layers.LayerSystem) – stratified medium
  • k_parallel (numpy.ndarray or str) – in-plane wavenumber. If ‘default’, use smuthi.fields.default_Sommerfeld_k_parallel_array
  • resolution (float or None) – spatial resolution of the lookup in the radial direction
  • cuda_blocksize (int) – threads per block when calling CUDA kernel
class smuthi.linearsystem.linear_system.CouplingMatrixVolumeLookup(vacuum_wavelength, particle_list, layer_system, k_parallel='default', resolution=None)

Base class for 3D lookup based coupling matrix either on CPU or on GPU (CUDA).

Parameters:
  • vacuum_wavelength (float) – vacuum wavelength in length units
  • particle_list (list) – list of sumthi.particles.Particle objects
  • layer_system (smuthi.layers.LayerSystem) – stratified medium
  • k_parallel (numpy.ndarray or str) – in-plane wavenumber. If ‘default’, use smuthi.fields.default_Sommerfeld_k_parallel_array
  • resolution (float or None) – spatial resolution of the lookup in the radial direction
class smuthi.linearsystem.linear_system.CouplingMatrixVolumeLookupCPU(vacuum_wavelength, particle_list, layer_system, k_parallel='default', resolution=None, interpolator_kind='cubic')

Class for 3D lookup based coupling matrix running on CPU. This is used when no suitable GPU device is detected or when PyCuda is not installed.

Parameters:
  • vacuum_wavelength (float) – vacuum wavelength in length units
  • particle_list (list) – list of sumthi.particles.Particle objects
  • layer_system (smuthi.layers.LayerSystem) – stratified medium
  • k_parallel (numpy.ndarray or str) – in-plane wavenumber. If ‘default’, use smuthi.fields.default_Sommerfeld_k_parallel_array
  • resolution (float or None) – spatial resolution of the lookup in the radial direction
  • interpolator_kind (str) – ‘linear’ or ‘cubic’ interpolation
class smuthi.linearsystem.linear_system.CouplingMatrixVolumeLookupCUDA(vacuum_wavelength, particle_list, layer_system, k_parallel='default', resolution=None, cuda_blocksize=None, interpolator_kind='linear')

Class for 3D lookup based coupling matrix running on GPU.

Parameters:
  • vacuum_wavelength (float) – vacuum wavelength in length units
  • particle_list (list) – list of sumthi.particles.Particle objects
  • layer_system (smuthi.layers.LayerSystem) – stratified medium
  • k_parallel (numpy.ndarray or str) – in-plane wavenumber. If ‘default’, use smuthi.fields.default_Sommerfeld_k_parallel_array
  • resolution (float or None) – spatial resolution of the lookup in the radial direction
  • cuda_blocksize (int) – threads per block for cuda call
  • interpolator_kind (str) – ‘linear’ (default) or ‘cubic’ interpolation
class smuthi.linearsystem.linear_system.LinearSystem(particle_list, initial_field, layer_system, k_parallel='default', solver_type='LU', solver_tolerance=0.0001, store_coupling_matrix=True, coupling_matrix_lookup_resolution=None, interpolator_kind='cubic', cuda_blocksize=None, periodicity=None, ewald_sum_separation_parameter='default', number_of_threads_periodic='default', use_pvwf_coupling=False, pvwf_coupling_k_parallel=None)

Manage the assembly and solution of the linear system of equations.

Parameters:
  • particle_list (list) – List of smuthi.particles.Particle objects
  • initial_field (smuthi.initial_field.InitialField) – Initial field object
  • layer_system (smuthi.layers.LayerSystem) – Stratified medium
  • k_parallel (numpy.ndarray or str) – in-plane wavenumber. If ‘default’, use smuthi.fields.default_Sommerfeld_k_parallel_array
  • solver_type (str) – What solver to use? Options: ‘LU’ for LU factorization, ‘gmres’ for GMRES iterative solver
  • store_coupling_matrix (bool) – If True (default), the coupling matrix is stored. Otherwise it is recomputed on the fly during each iteration of the solver.
  • coupling_matrix_lookup_resolution (float or None) – If type float, compute particle coupling by interpolation of a lookup table with that spacial resolution. A smaller number implies higher accuracy and memory footprint. If None (default), don’t use a lookup table but compute the coupling directly. This is more suitable for a small particle number.
  • interpolator_kind (str) – interpolation order to be used, e.g. ‘linear’ or ‘cubic’. This argument is ignored if coupling_matrix_lookup_resolution is None. In general, cubic interpolation is more accurate but a bit slower than linear.
  • periodicity (tuple) – tuple (a1, a2) specifying two 3-dimensional lattice vectors in Carthesian coordinates with a1, a2 (numpy.ndarrays)
  • ewald_sum_separation_parameter (float) – Used to separate the real and reciprocal lattice sums to evaluate particle coupling in periodic lattices.
  • number_of_threads_periodic (int or str) – sets the number of threats used in a simulation with periodic particle arrangements if ‘default’, all available CPU cores are used if negative, all but number_of_threads_periodic are used
  • use_pvwf_coupling (bool) – If set to True, plane wave coupling is used to calculate the direct. Currently only possible in combination with direct solver strategy.
  • pvwf_coupling_k_parallel (array) – k-parallel for PVWF coupling
compute_coupling_matrix()

Initialize coupling matrix object.

compute_initial_field_coefficients()

Evaluate initial field coefficients.

compute_t_matrix()

Initialize T-matrix object.

prepare()
solve()

Compute scattered field coefficients and store them in the particles’ spherical wave expansion objects.

class smuthi.linearsystem.linear_system.MasterMatrix(t_matrix, coupling_matrix)

Represent the master matrix \(M = 1 - TW\) as a linear operator.

Parameters:
class smuthi.linearsystem.linear_system.SystemMatrix(particle_list)

A system matrix is an abstract linear operator that operates on a system coefficient vector, i.e. a vector \(c = c_{\tau,l,m}^i\), where \((\tau, l, m)\) are the multipole indices and \(i\) indicates the particle number. In other words, if we have a spherical wave expansion for each particle, and write all the expansion coefficients of these expansions into one (long) array, what we get is a system vector.

index(i, tau, l, m)
Parameters:
  • i (int) – particle number
  • tau (int) – spherical polarization index
  • l (int) – multipole degree
  • m (int) – multipole order
Returns:

Position in a system vector that corresponds to the \((\tau, l, m)\) coefficient of the i-th particle.
index_block(i)
Parameters:i (int) – number of particle
Returns:indices that correspond to the coefficients for that particle
class smuthi.linearsystem.linear_system.TMatrix(particle_list)

Collect the particle T-matrices in a global lienear operator.

Parameters:particle_list (list) – List of smuthi.particles.Particle objects containing a t_matrix attribute.
right_hand_side()

The right hand side of the linear system is given by \(\sum_{\tau l m} T^i_{\tau l m} a^i_{\tau l m }\)

Returns:right hand side as a complex numpy.ndarray

linearsystem.linear_system_cuda

This module contains CUDA source code for the evaluation of the coupling matrix from lookups.